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ABSTRACT 
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S-) \ Cosmic ray nuclei fluxes are expected to be measured with high precision in the 

near future. For instance, high quality data on the antiproton component could give 
important clues about the nature of the astronomical dark matter. A very good un- 
derstanding of the different aspects of cosmic ray propagation is therefore necessary. 
In this paper, we use cosmic ray nuclei data to give constraints on the diffusion pa- 
rameters. Propagation is studied with semi-analytical solutions of a diffusion model, 
and we give new analytical solutions for radioactively produced species. Our model 
includes convection and reacceleration as well as the standard energy losses. We 
perform a \ 2 analysis over B/C data for a large number of configurations obtained 
by varying the relevant parameters of the diffusion model. A very good agreement 
with B/C data arises for a number of configurations, all of which are compatible with 
sub-Fe/Fe data. Different source spectra Q(E) and diffusion coefficients K{E) have 
been tried, but for both parameters only one form gives a good fit. Another impor- 
tant result is that models without convection or without reacceleration are excluded. 
We find that the various parameters, i.e. the diffusion coefficient normalisation Kq 
and spectral index 5, the halo thickness L, the Alfven velocity V a , and the convection 
velocity V c are strongly correlated. We obtain limits on the spectral index S of the 
diffusion coefficient, and in particular we exclude a Kolmogorov spectrum (<5 = 1/3). 
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1. Introduction 

Understanding the composition and spectral features of cosmic rays has always been an 
astrophysical challenge. On one hand, the observational data have long been scarce and suffered 
from large uncertainties. On the other hand, the theoretical predictions to which these data 
should be compared to have also suffered from several drawbacks. Composition and spectra 
arise from the nuclear interaction of an initial distribution of energetic particles with interstellar 
matter (spallations) and their electromagnetic interactions with galactic magnetic fields (acceler- 
ation and diffusive reacceleration) . First, the nuclear cross sections to be used were not very well 
known until recently. Second, our knowledge of the galactic magnetic field is far from complete. 
Cosmic rays are sensitive to its scale inhomogeneities (diffusion) which are not well observed. 
They are also sensitive to the presence of plasma shock-waves (acceleration in localized sources 
and diffusive reacceleration). Third, composition and spectra are altered as the cosmic rays 
enter the solar magnetic field, so that some more modelling has to be done in order to infer 
interstellar spectra from observations. 

However, despite these sources of uncertainty, some gross features of the cosmic ray prop- 
erties are well established. First, the study of secondary-to-primary ratio shows that a cosmic 
ray has crossed an average of about x = 9 g cm -2 of interstellar matter between its acceleration 
and its detection. Second, the isotopic ratio of radioactive species shows that it took about 
r = 20 Myr between the same events. As cosmic rays have a velocity close to c, we can infer the 
average density of the medium to be n = x/(t cmn) ~ 0.3 cm -3 . As the disk density is about 
1 cm~ 3 , the cosmic rays must spend a fraction of time in an empty region, called the diffusion 
halo. Cosmic rays are produced and destroyed in the galactic disk, they diffuse in a larger zone 
whose characteristics are not well known, and eventually they can escape from this zone. To 
have a good modelling of cosmic ray propagation, we should know the geometry and size of the 
diffusion zone, the characteristics of the galactic magnetic field, and the sources. 

In this paper, we use the existing data on nuclei, gathered by balloon-borne and space 
experiments during the last thirty years, to put constraints on the parameters describing the 
propagation of cosmic ray nuclei. One consequence is that many of the large uncertainties 
affecting the calculation of the antiproton flux could be strongly reduced (F. Donato & al., in 
preparation (12)). This is of utmost relevance for the study of exotic sources of antiprotons and 
antideuterons. 

2. The diffusion model 

It has been recognized for a long time that the relevant physical propagation model to be 
used is the diffusion model (Berezinskii & al. 1990; Ginzburg &: Syrovatskii 1964), though the 
so-called leaky box model has been widely preferred for decades because of its simplicity. 

The steady-state differential density N 3 (E, r) of the nucleus j as a function of energy E 
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and position r in the Galaxy, is given by (see for instance Berezinskii & al. 1990) 

V. { K>VW - V c N j ) - A (Yf Ek + IL { » NS) _ 1 + =qi+ -£ f«iV* 

v 1 k / / m. k >mj 

The first terms represent diffusion (iP is the diffusion coefficient) and convection (V c is the 
convection velocity). The divergence of this velocity, expressed in the next term gives rise to an 
energy loss term connected with the adiabatic expansion of cosmic rays. Further, we have to take 
into account ionization and coulombian losses (the energy loss coefficients are specified below), 
plus a reacceleration term in first order derivative (all included in b 3 ) and finally a second order 
derivative in E for the associated second order term in reacceleration (d 3 is the energy diffusion 
coefficient). These stand for the continuous losses. The last term of the l.h.s. takes care of the 
disappearance of the nucleus j (T 3 for short) due to its collisions with interstellar matter (ism). 
In the r.h.s., the source term q 3 takes into account the primary production and acceleration of 
nuclei described by an injection spectrum (for the sake of clarity, we have not written down 
the terms describing the contribution of radioactive species). Finally, the last term is for the 
secondary j sources, namely spallation contribution 2 T k3 from all other heavier nuclei. We use 
here the notation E for total energy and E^ for kinetic energy. 

One needs to solve a complete triangular-like set of coupled equations since only heaviest 
nuclei contribute to a given nucleus. Quantities in this equation are functions of spatial coordi- 
nates (not time, steady-state being assumed) and of energy. Unless stated otherwise, the word 
energy will stand for kinetic energy per nucleon (T), since this is the appropriate parameter to 
be used, as it is conserved in spallation reactions (see footnote 2). 

The other very popular model is the leaky box, in which all the quantities are spatially 
averaged (so that convection has no meaning), and the diffusion term is then replaced by an 
escape term 

-V.(K 3 VN 3 ) < — > N 3 /r esc 

which has the meaning of a residence time r esc (Myr) in the confinement volume. Note that 
in this case, the mean time (r) spent in the box and the mean column density (x) (in units 
of g cm -2 ) of crossed matter are treated on the same footing since the model is homogeneous. 
This is not true for multi-zone models and especially for the diffusion model where cosmic rays 
spend most of their time in the diffusion halo which has a zero matter density. It may seem very 
surprising that this simplified model is able to reproduce many observations on stable nuclei. 
Actually, it can be shown that leaky box models are often equivalent to diffusion models, as far 
as stable species are concerned (for a generic discussion on the "leakage-lifetime" approximation 
in cosmic ray diffusion, see Jones 70). This is not obvious from the equations, but it can be 
seen more readily if we solve them in a formalism called the weighted slab technique. Basically, 



2 We use here the so-called straight-ahead approximation which implies the conservation of kinetic energy per 
nucleon during a spallation process 

p oo p oo 

/ n H v'N k {T')a k3 (T,T')dT' = n„v' N k (T')a kj (T)S(T - T')dT' = n H vN k (T)a kj (T) (2) 
Jo Jo 
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it consists in writing a general solution under the form (see for example Berezinskii et al. 1990; 
Ginzburg & Syrovatskii 1964) 

POO 

N j (r) = / N j {x)G(r,x)dx (3) 
Jo 

Thus, (1) decouples into two independent equations. The first one involves G(r,x). It depends 
only on the geometry of the problem and on the chosen diffusion scheme (or the leaky box 
parameters), but not on the particular species j. The other one involves N 3 (x) and contains the 
nuclear physics aspects of the propagation. Different models (diffusion or leaky box) correspond 
to different G(r,x), but the equations on N J (x) are the same. G is the path length distribution 
and is interpreted as the distribution of probability that a nucleus j goes through a column 
density x before reaching Earth. The global solution is then obtained by a convolution of the 
surviving quantity N^(x) under nuclear processes after x g cm -2 , with the weighted probability 
associated. 

In a leaky box model, it can be shown that the normalized path length distribution function 
is given by 

G(x) = lexp(f) (4) 

^esc \ ^csc / 

where the average quantity of matter crossed by a cosmic ray is given by 

/•oo 

(x) = / xG(x) dx = A esc (5) 
Jo 

In the diffusion model, it has been shown (Owens 1976) that for a wide class of geometries, 
the function G is given by an infinite series of exponentials involving the diffusion coefficient. 
Moreover, the physical parameters are such that for sufficiently large grammage x, only the first 
exponential contributes so that we recover a leaky box model (Jones 1970; Ginzburg, Khazan 
& Ptuskin 1980; Schlickeiser & Lerche 1985; Lerche & Schlickeiser 1985). For the model used 
in this paper, it can be shown that A esc = fhnvhL / K in the limit L <C R. However, this is no 
longer true for lower x, and then the leaky box model fails. This may explain, in conjunction 
with poorly known cross sections for heavy nuclei, why for many years cosmic ray physicists 
had to use modifications (Lezniak & Webber 1979; Garcia-Munoz Sz al. 1987; Webber & al. 
1998b) of the leaky box path length distribution at small grammage to reconcile observations 
of B/C with those of sub-Fe/Fe. Even worst, this approximate equivalence between leaky box 
and diffusion model breaks down for radioactive species (Prishchep & Ptuskin 1975). 

This motivates us to choose the diffusion model and solve the corresponding equation (1). 
From a more general point of view, note also that the weighted slab technique mentionned above 
should be used carefully. Actually, it was recently realized (Ptuskin, Jones & Ormes 1996) that 
pre-1996 estimations of cosmic ray parameters with this technique were not exact. First, at 
fixed energy per nucleon, the rigidity and therefore the propagation features depend on the 
nuclear species, whatever the energy. This prevents the above-mentionned decoupling. Second, 
energy losses that strongly depend on A and Z also take place at low energy. Ptuskin et al. 
1996 presented a modification of this formalism which is correct in both the highly and non- 
relativistic regimes. These fundamental findings were confirmed by a more complete numerical 
treatment (Stephens & Streitmatter 1998) and give rise to differences of as much as 20%. 
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For the nuclear part, we adopted the direct approach of the shower technique. This means 
that the flux is first evaluated for the heavier primary cosmic ray, for which the diffusion equation 
is simple and does not couple to any other species. Then, the flux of the next nucleus is 
computed, with a spallation term depending only on the heavier nucleus, whose flux is known 
from the previous step. This procedure is repeated for all the nuclei down to the lightest one. 
We started at Z = 30, the heavier non negligeable primary found in cosmic radiation (Binns & 
al. 1989), since (Z > 30) < 10 -3 Fe. For all light nuclei, we checked that it is sufficient to start 
from Z = 16 (S). For the detail of solutions, the reader is referred to Appendix A. The next 
section is devoted to the inputs of the model. 



3. The parameters of the model 

3.1. Geometry of the Galaxy 

The geometry of the problem used here is a classical cylindrical box (see for example Webber 
et al. 1992) whose radial extension is R, with a disk of thickness 2h and a halo of half-height L. 
Only L is not fixed (see §5) and R = 20 kpc, h = 100 pc. Sources and interactions with matter 
are confined to the thin disk and diffusion which occurs throughout disc and halo with the same 
strength is independent of space coordinates. The solar system is located in the galactic disc 
[z = 0) and at a centrogalactic distance Rq = 8 kpc (Stanek & Garnavich 1998; Alves 2000). 



3.2. Cross sections 



Cross sections play a crucial role in propagating the cosmic rays throughout the Galaxy. We 
have to distinguish bewtween the total inelastic cross section and the spallation cross section. 



3.2.1. Total inelastic cross section 



The total inelastic cross section, which is actually a reaction cross section, is defined by 

—tot ~tot ~tot (c\ 

a inel — a ~ ° elastic \°J 

Various empirical modifications of the original (Bradt & Peters 1950) geometrical approach 
have been used to parameterize this total cross section. Letaw, Silberberg & Tsao (1983) 
and Silberberg & Tsao (1990) produced the basic energy dependent equations further refined 
in Sihver & al. (1993) and in Wellish & Axen (1996). Recently Tripathi, Cucinotta & Wilson 
(1997a, 1997b, 1999) proposed a universal parameterization which is valid for all nucleus-nucleus 
reactions, including neutron induced reaction. All of these approaches have been compared in 
Silberberg, Tsao & Barghouty (1998), where the preference were given to Tripathi et al. (1997a, 
1997b, 1999). For this reason, we adopted the latter parameterization throughout this study. 
The precision can be estimated to be about 5%. 
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3.2.2. Spallation cross section 

The fragmentation cross sections are somewhat more complicated. Basically, we can sepa- 
rate three distinct approaches : (i) microscopical description (Bondorf & al. 1995; Ramsey & 
al. 1998) (ii) semi-empirical formulae (Silberberg et al. 1998; Tsao, Silberberg & Barghouty, 
1998; Tsao et al. 1999) and (iii) empirical formulae (Webber, Kish & Schrier 1990d). The first 
one takes into account the fundamental nuclear physics but is very computer time consuming 
and has not reached a precision as good as the two other approaches. The second method 
takes into account more phenomenological inputs, and is better designed for the unmeasured 
fragmentation cross sections which correspond mainly to the ultra heavy {Z > 30) nuclei. The 
third one breaks down in this region, but can be used for 3 < Z < 30 because it is adapted to fit 
measurements. Consequently, we used the code of Webber et al. (1990d) available on the web 3 
updated with new parameters (table V of Webber & al. 1998a). This parameterization takes 
advantage of the fact that about 98% of all the reaction involved in cosmic ray propagation have 
now at least one point of energy measured. This set of formulae gives a precision of about 10% 
but can not be further refined (Webber et al. 1998a). The code is extended for spallation on He 
with the parameterization given in Ferrando et al. (1988) 4 . 

To go further than this treatment, one should keep in mind the following points. First, 
the process of single or double nucleon removal is an important channel for the production 
of numerous nuclei. It has been parameterized in detail in Norbury & Townsend (1993) and 
in Norbury & Mueller (1994), although in the following we will consider Webber's approach 
sufficient for our purpose. Second, most of the parameterizations make use of the straight- 
ahead approximation (see footnote 2). Tsao & al. (1995) relaxed this condition and find a 5% 
effect around the 1 GeV peak of B/C (note that this effect is negligable for primary to primary 
and secondary to secondary ratios). Finally, it was found that all the measurements are not 
always consistent and systematic effects are likely to be present (Taddeucci & al. 1997; Vonach 
& al. 1997; Zeitlin & al. 1997; Flesh & al. 1999; Korejwo & al. 1999). 

To summarize, we consider that Webber's code is very well suited for modelling spallation 
cross sections, thanks to the many experiments developed these last years (Webber, Kish &: 
Schrier 1990a, 1990b, 1990c; Webber et al. 1998b, 1998c; Chen et al. 1997). The overall 
precision can be roughly estimated to be better than 10% on H, and better than 20% on He. 

3.3. Nuclear parameters 

Apart from the stable species, we can distinguish two categories of nuclei to be propagated, 
depending on whether they live long enough to be observed or they decay quickly into other 
nuclei. 



3 http://spdsch.phys.lsu.edu/SPDSCHJ D ages/Software_Pages/Cross_Section_Calcs/CrossSectionCalcs.html 

4 The parameterization given in this paper should be modified at low energy and for large AZ (P. Ferrando 
2001, private communication). 
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3.3.1. Ghost nuclei 

Ghost nuclei are defined as those having a half-lifetime larger than ~ 1 ms (so that they 
are detected in cross section measurements) and less than few hundred years (so that we can 
consider that disintegration is fast on a propagation timescale). Several reaction channels may 
yield a given nucleus. This means that for this nucleus, we have to add to the direct production 
the contributions of all the short lived parent nuclei 

ghost 

effective = ^ . + £ ak ^ xB r(X^j) (7) 

nuclei X 

where Br(X — > j) is the branching ratio of the corresponding channel 5 . A table of all the 
disintegration chains can be found in Letaw & al. 1984. To see the effect of the recent discovery 
of new intermediate short-lived nuclei and of the changes in branching ratio, we reconstructed 
from a recent compilation of nuclear properties (Audi et al. 1997) all the ghost nuclei and 
compared them with the Letaw table. We found that the net effect is negligible. This is because 
the discovered radioelements are very far from the valley of stability and their production is 
suppressed (as it can be observed experimentally). Such reaction chains should accordingly be 
taken as a fixed input for the future. 



3.3.2. Radioactive nuclei 

Radioactive nuclei can be separated into three groups : unstable under (3 decay, unstable 
under electronic capture and unstable under both reactions. All the physical effects related to 
these nuclei will be analysed in detail in a separate paper (F. Donato & al., in preparation (13)). 



3.4. Diffusion coefficient 

Physically, diffusion arises because charged particles interact with the galactic magnetic 
field inhomogeneities. This is an energy-dependent process because higher energy particles are 
sensitive to larger spatial scales. As a result, the diffusion coefficient is related to the power 
spectrum of these inhomogeneities, which is poorly known. Several analytical forms of this 
energy dependence have been assumed in the literature. In particular, leaky box models give a 
purely phenomenological form 

A csc = A /3 if n < 1Z C in GV 

A esc = X P {K/K c r 5 if K > K c in GV 



5 To give an example, such a reaction reads for %Be 



<"lB7 e = °i^lBe + 49.2% a k ^ Lt + 4.1% a k ^ Li 



(8) 
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where 1Z = p/Z stands for the particle rigidity. A more fundamental origin may be found for 
the previous expression, as one can show in a one-dimensional diffusion model with convection 
(Jones & al. 2001). All the effects that led to that phenomenological form are explicitely taken 
into account in our treatment and we may therefore consider a more natural form for K(E), 
inferred from magnetohydrodynamics considerations (Ptuskin et al. 1997): 

K(E) = K (3x1l 5 (9) 

where the normalisation Kq is expressed in kpc 2 Myr -1 . Nevertheless, we have also tested leaky 
box inspired forms, but these never give good fits to the data (see next section). 

3.5. Sources 

A charged particle is called a cosmic ray when it has been accelerated to an energy greater 
than 100 MeV/n as far as our study is concerned. Several physical processes may be responsible 
for this energy gain; they all involve an interaction of the particle with a plasma shockwave front. 
This is responsible for the power-law energy dependence (see below). The term source refers to 
the place where the charged particles have been promoted to cosmic rays. We present here the 
various hypotheses about their spatial and spectral distributions and their compositions. 

3.5.1. Spatial distribution 

First, as explained above, we assumed that the Galaxy has a cylindrical symmetry and that 
the sources are located in the galactic disc. So we only have to specify a radial distribution. We 
have used the Case & Bhattacharya (1996, 1998) distributions, but the results presented below 
are not sensitive to this particular choice. 

3.5.2. Spectral distribution 
The source spectral distributions can be split in two terms 

q\E)=qiQ3(E) (10) 

where represent the composition and Q 3 {E) the energy dependence. Several spectra have 
been used in the past (Engelmann & al. 1985), 

Q j (E) oc p~ a i 

Qi(E) oc PE~% (11) 
Q\E) oc E~? 

where p is the momentum and where the spectral index ctj may depend on the species j consid- 
ered. These forms are equivalent at high energies but we focused on the first one only (equivalent 
to a power-law in rigidity) as the others do not yield good fits to the data. The sum 5 + ay 
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gives the overall power-law index of each species at high energy (~ 1 TeV/n): this is because 
the measured flux at these energy is proportional to Qi(E)/K(E) oc lZ &+a i . In other words, 
at sufficiently high energy, energetic redistribution and spallation reactions do not affect prop- 
agation. Consequently, measured spectral indexes enter as a fixed input, and we took them in 
Wiebel-Sooth, Biermann & Meyer (1998). For consistency, we required that whatever 5, the 
parameters ay are adjusted so as to reproduce the above-mentionned values. This ensures that 
the spectra we compute agree with data at high energy. 



and data. Alternatively, a more fundamental approach can be used. It may be taken either as 
the solar system composition convoluted with the fip (First Ionisation Potential), or convoluted 
with the volatility. As pointed out in Webber (1997), the two alternatives are not so different. 
In this work, the isotopic composition of each element was taken from Anders & Grevesse (1989) 
and Grevesse & Sauval (1998), and the q$ dependence of the fip has been taken from Binns et 
al. (1989). Then the primary cosmic ray elemental composition has been adjusted to fit the 
heao-3 data at 10.6 GeV/n, while keeping constant the relative isotopic abundances. 



Once we have chosen a generic model, we have to include effects which affect the spectra at 
low energies. A first step is the inclusion of energetic losses. Then a most complete treatment 
must take into account reacceleration and convection. 



There are two types of energy losses which are relevant for nuclei: ionization losses in the 
ISM neutral matter (90% H and 10% He), and Coulomb energy losses in a completely ionized 
plasma, dominated by scattering off the thermal electrons. We will use < n e >~ 0.033 cm -3 , 
and T e ~ 10 4 K (Nordgren 1992). Complete formulae for these two effects are compiled, for 
example, in Strong & Moskalenko (1998) or in Mannheim & Schlickeiser (1994). There is also a 
third type of loss related to the existence of a galactic wind. We recall that in equation (1), there 
is a term proportional to V.14. As we choose a constant wind velocity in the z direction (see 
farther), we could conclude at a first glance that this term vanishes. Actually, this is the case 
everywhere in the Galaxy except at z = where a discontinuity occurs, due to the opposite sign 
of the wind velocity above and below the galactic plane. One gets a term that can be expressed 
in the same form as ionization and coulombian losses with an effective term 



3. 5. 3. Composition 



The composition q$ is usually inferred from a comparison between leaky box calculations 



3.6. Galactic wind, energy loss, reacceleration 



3.6.1. Energy losses 




(12) 
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Ek stands for the total kinetic energy and it should not be confused with the kinetic energy per 
nucleus frequently used in this paper. Finally, because of these effects the propagation equation 
must be solved numerically. 



3. 6. 2. Reacceleration 



As the fractional energy changes in a single collision are small, we can treat Fermi accel- 
eration (Fermi 1949; 1954) using the Fokker-Planck formalism (Blandford & Eichler 1987). In 
one dimension 



dt dx 



d_ 

dp 



1 d 



At 1 1 2 dp 



^(Ap) 2 



At 



f 



(13) 



where (Ap/At) and ((Ap) 2 / At) are expressed in terms of the probability for changing the 
momentum p to p + Ap during time At. A simplification occurs when the recoil of the scatterer 
can be ignored. In general, the principle of detailed balance ensures that the probability for a 
gain is equal to that for a loss. This implies 



Ap 
At 



which reduces the equation into 



dt dx 



where 



K. 



pp 



1 d I (Ap) 2 
2dp \ At 



1*4 (/) 



1 /{Ap) 2 

2 \ At 



(14) 



(15) 



(16) 



From another point of view, we can start with a collisionless Boltzmann equation and expand 
it up to second order in perturbed quantities (magnetic irregularities). It gives the full diffusion 
equation in quasi-linear theory (Schlickeiser 1986; Kulsrud & Pearce 1969). Assuming a constant 
flow of magnetic irregularities (corresponding to a constant galactic wind) in the z direction, we 
find 

2 df] 



f = -y c V/ + V(*V/ ) + ^| 



K, 



ppP dp 



+ Q(r,p,t) 



(17) 



which is the propagation equation we used in Appendix A (up to several nuclear terms, cf 
eq. (Al)). In this theory, the momentum diffusion coefficient K pp is naturally connected to the 
space diffusion coefficient K by the relation 



K, 



pp 



V 2 
V 

m 



(18) 



where V a is the Alfven speed of the scattering centers. 



We emphasize that expressions (15) and (17) are just simplifications of the full Fokker- 
Planck equations. As pointed out by Ostrowski & Siemieniec-Oziehlo (1997), if this simplification 
is natural in the kinetic description of a gas of scattering particles, this is no longer the case 
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for particles scattered off external 'heavy' scattering centers. Nevertheless, this approximation 
which arises in a number of astrophysical cases is sufficient for our purpose. We thus have 
to treat a simple energy diffusion coefficient, which we evaluated in the no-recoil hard sphere 
scattering centers approximation as 

where o~K n K describes the rate of collisions with scatterers which is related to the mean free 
path by Xk = I/o-r^k- As K{E) = (1/3)\kv, the previous relation gives the final coefficient 
to be used in our calculation 

K »-\ v -f$i m 

Our approch is basically similar to that of Heinbach & Simon (1995). The meaning of the 
Alfvenic speed derived in their analysis is not very clear, since the numerical factor in front 
of expression (20) may vary according to the hypotheses made for the scattering process. The 
assumption that reacceleration occurs only in the thin disc is supported by recent complete 
magnetohydrodynamics simulations (Ptuskin et al. 1997) describing the evolution of the hot 
gas, cosmic rays and magnetic field. 



3.6.3. Integration into equation 

We define 

7 / x /dE\ I dE\ /dE\ 

\ aT I Ion \ al I Coul \ al I Adiab 

and with K J EE (E) given by (20), one obtains (see Appendix) 

A{Ni(0) = & - 2h^L j&L s (*W(0) - Ki E (E)-^Ni(0)} (22) 

This is a second order equation which is solved for each order i of the Bessel decomposition. 
A\ is given by (A10) and Q? (A5) is a generic source term (see Appendix A for details). Rela- 
tion (22) has been developed on our energy array so as to be transformed into a matrix equation 
that is directly inverted. Two boundary conditions need to be implemented. The effects of 
energy losses and diffusive reacceleration have been assumed to be negligible at the high en- 
ergy tip (100 GeV/n) of the spectrum. We have furthermore imposed a vanishing curvature 
(d 2 Nl /dE 2 = 0) at the lowest energy point. 



3.6.4- Galactic wind 

It has been recognized for a long time that a thin disc configuration would be disrupted by 
cosmic ray pressure (Parker 1965, 1966; Ko, Dougherty &; McKenzie 1991). It can be stabilized 
by the presence of a halo, but further considerations imply that this halo would not be static 
either (presence of a convective wind). Consequences of a wind has been firstly investigated 
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by Ipavitch (1975) and since then it has been observed in other galaxies (Lerche & Schlickeiser 
1982a, 1982b; Reich & Reich 1988; Werner 1988). In our own, its effects on nuclei have been 
investigated in various models (Jones 1979; Kota & Owens 1980; Freedman & al. 1980; Bloemen 
& al. 1993), and in this work we adopted a very simple and tractable form for V c throughout 
the diffusive volume (dV c /dz = 0), following Webber, Lee & Gupta (1992). A constant wind 
with free escape boundary is a reasonable approximation to a more sophisticated variable-wind 
model (see also conclusions from magnetohydrodynamics simulations (Ptuskin et al. 1997)). 

3.6.5. Comparison with other works 

Our semi-analytical two-dimensional model lies, in some sense, between the complete nu- 
merical resolution of Strong & Moskalenko (1998) and the recent study of Jones et al. (2001). 
As compared to the first one, we have a simplified description of the matter distribution, for 
which the full analysis of the parameter space does not require too much CPU time. The second 
approach can be considered as the limiting case for our two-dimension model when the radial 
extension is much larger than the halo size (R^> L). Note that a similar analysis was performed 
by Letaw, Silberberg &: Tsao (1993) in the framework of a leaky box model or Seo & Ptuskin 
(1994) in a one-dimensional diffusion model but their less systematic approach does not allow 
a full exploration of parameter space. 

4. Data Analysis 

For the aim of our analysis, we can consider different classes of flux ratios: primary-to- 
primary (e.g. C/O), secondary-to-primary (e.g. B/C or sub-Fe/Fe), secondary-to-secondary 
(e.g. Li/B or Be/B), ratio of isotopes either stable (e.g. 10 B/ 11 B) or unstable (e.g. 10 Be/ 9 Be). 
Each of these may be an indicator of some dominant physical phenomenon and be particularly 
sensitive to the corresponding diffusion parameters. The ratio of two primaries is practically 
insensitive to changes in all the parameters, since they have the same origin and undergo the 
same physical processes. So, it can be a very useful tool to fix their source abundances. We will 
return to this subject. Similar conclusions, even if less strong, may be drawn for the ratio of two 
isotopes of the same species, such as 10 B/ 11 B. Indeed, at very low energy values this quantity is 
slightly affected by changes in the injection spectra, but the effect is too weak to constrain free 
parameters. 

The most sensitive quantity is B/C, as B is purely secondary and its main progenitor C 
is primary. The shape of this ratio is seriously modified by changes in the diffusion coefficient, 
in the height of the diffusion halo, in the convective and alfvenic velocities. Moreover, it is 
also the quantity measured with the best accuracy, so that it is ideal to test models. Indeed, 
being the ratio of two nuclei having similar Z, it is less sensitive to systematic errors than single 
fluxes or other ratios of nuclei with more distant Z. For the same reasons, the sub-Fe/Fe may 
also be useful. Unfortunately, since existing data are still affected by sizeable experimental 
errors, we can only use them to cross-check the validity of B/C but not to further constrain 
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the parameters under scrutiny. Another particularly interesting quantity is the ratio 10 Be/ 9 Be. 
Since 10 Be is a radioactive element, it is very sensitive to the processes which can occur in the 
halo. Therefore, we could use this ratio in particular to bracket the size of the halo (F. Donato 
& al., in preparation). 



4.1. Solar modulation 

We have estimated the effect of the solar wind on the energies and intensities of cosmic-rays 
following the prescriptions of the force field approximation (Perko 1987). The modification in 
the total interstellar energies of a nucleus with charge Z and atomic number A, corresponds to 
the shift Zcj) (we set the electronic charge equal to unity): 

E TOA /A = E ls /A - \Z\<P/A. (23) 

Here E TOA and E ls correspond to the top-of-atmosphere (modulated) and interstellar total 
energy, respectively. The solar modulation parameter <f> has the dimension of a rigidity (or an 
electric potential), and its value varies according the 11-years solar cycle, being greater for a 
period of maximal solar activity. Often people refer to the equivalent quantity $ = \Z\ <j>/A — 
i^>. Once the momenta at the Earth p TOA and at the boundaries of the heliosphere p ls are 
determined, the interstellar flux of the considered nucleus is related to the TOA flux according 
to the simple rule: 

* T ° A (£ T ° A ) (P^\ 2 (24) 
3>is (e 1s ) \ pis J ^ 

The determination of the modulation parameter <1> is totally phenomenological and suffers from 
some uncertainty. As explained in the following, we will deal with data taken around period of 
minimal solar activity, for which we fixed <I> = 250 MV. 



4.2. The dataset we used 

In order to test our diffusion model we have widely employed the data taken by the exper- 
iment which was launched on board the NASA HEAO-3 satellite (Engelmann et al. 1990). The 
relative abundances of elements, with charge from 4 to 28 and for energy ranging from 0.6 to 
35 GeV/n, have been measured with unprecedent accuracy. Data have been taken in 1979-80, 
around a solar minimal activity. For the case of B/C, heao-3 quoted errors are 2-3%. To 
perform the analysis, we also included data from balloons (Dwyer & Meyer 1987) and from the 
ISEE-3 experiment (Krombel & Wiedenbeck 1988), even if the relevant error bars are wider. 
The first one collected data in 1973-75 for energies spanning from around 1.7 to 7 GeV/n. The 
second experiment was operating during 1979-81 on board a spacecraft, in the energy range 
100-200 MeV/n. In some of the figures presented below, we also plot - for purely illustrative 
goals - the data point from IMP-8 (Garcia-Munoz et al. 1987) and the voyager experiments 
(Lukasiak, McDonald & Webber 1999) (we did not include Ulysse data point (DuVernois & al. 
1996), since it corresponds to a period of maximal solar activity). They were not included in 



-14- 



the x 2 analysis since they are the most sensitive to solar modulation. Nevertheless, we have 
checked that the best x 2 values were not significantly modified when these points were added. 

As regards the sub-Fe(Sc+V+Ti)/Fe ratio, we used data from HEAO-3 (Engelmann et al. 
1990) and from balloons (Dwyer & Meyer 1987). In both cases the error bars, around 10%, are 
significantly larger than for B/C. 



We varied the relevant parameters Kq, L, V c , V a , and 5 of the diffusion model described 
above. The constraints are much simpler to express with the combinations L, Kq/L, V c and 
V a j \JKq (this last expression appears naturally in Kee(E) describing reacceleration) . As for 
the primaries, we have adjusted the source abundance for nitrogen which is not a pure secondary. 
We calculated the x 2 over the 26 experimental points listed in 4.2, for each possible combination 
obtained varying the free parameters in the whole parameter space. In the following, we will 
focus on B/C and we cross-check the compatibility of the obtained parameters when applied to 
other fluxes. 



In a first step, the power 5 of rigidity in the expression for the diffusion coefficient (see 
Eq. (9)) has been fixed to 0.6. This will be the reference value for the bulk of our analysis. 
We only show models giving a x 2 l ess than 40. The best models have x 2 ~ 28. Some of these 
configurations have a very small halo size in which case the condition h -C L is no longer valid. 
Thus we also required L > 1 kpc in the whole analysis. 

The first result is that only the source spectrum q>{E) = qlp~ aj gives a good fit to the 



data. It will be used all the way long. For each L and Kq/L, we vary V^/V^o and V c and the 

best x 2 are shown in Fig. 1. All the good models are located in a narrow strip in the plane L 
- Kq/L, showing a strong correlation between these two parameters. However, there is still a 
degeneracy. It should be lifted with further specific analysis of radioactive nuclei (F. Donato 
&: al., in preparation). Conversely, for each V a /\/Ko and V c we varied L and Kq/L and the 
best values for the x 2 are depicted in Fig. 2. We see that they are gathered in a narrow range, 
namely: 8.5 < V c < 12 km sec" 1 and 410 < K/v 7 ^) < 530 km s" 1 / Mpc Myr" 1 / 2 . We do 
not find any model having a good x 2 without convection (V c = 0) or without reacceleration 
(V a = 0). We emphasize the fact that these four parameters are strongly correlated so that, 
when one of them is determined (for example L with radioactive nuclei), the allowed ranges for 
the others will be narrower than what could naively appear from the above two figures. 

The observed sub-Fe/Fe flux ratio does not give very powerful constraints. Nevertheless, all 
the models which give a good fit to B/C data are consistent with sub-Fe/Fe data. Figures 3 and 
4 illustrate these results with the particular choice L = 9.5 kpc, Kq/L = 0.00345 kpc Myr -1 , 
V c = 10.5 km s _1 and V a /\fK~Q~ = 470 km s^ 1 / Mpc Myr~ 1/2 , giving a reduced xl ~ 1-3. Note 
that all the configurations lying in the allowed ranges of Fig. 1 and 2 will provide very similar 
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Results 



5.1. Results of B/C and sub-Fe/Fe analysis for 5 = 0.6 
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curves. 

Up to now we have dealt with flux ratios, which partially hide the spectral features of the 
fluxes themselves. So we have also checked that the above configurations were consistent with 
the oxygen flux as measured by heao-3. In Fig. 5 the oxygen flux has been calculated for the 
same configuration as in Fig. 3 and 4. We can notice that the agreement with heao-3 data is 
not so satisfactory. The parameter which mainly prevents us from obtaining a good fit is the 
source spectrum. We recall that we set the sum 5 + ctj to the power-low index as determined 
in Wiebel-Sooth et al. (1998). Indeed, this number for the oxygen species is equal to 2.68 as 
derived from heao-3 and many other measurements. If instead we fix 5 + aoxygen to be 2.8 (as 
in Swordy & al. (1993)), the fit improves (see Fig. 5). 

5.2. Variation of the diffusion coefficient spectral index 

We also tested different values for 5, and we find that correspondingly the diffusion param- 
eters giving a good fit to B/C change. As an example, for a fixed value of the halo thickness 
L = 3 kpc, we find that 5 is allowed to vary between 0.5 and 0.84. This is displayed in Fig. 6. In 
the whole parameter space, the range of 5 extends from approximately 0.45 to 0.85, as one can 
see in Fig. 7. In particular the value 5 = 0.33 corresponding to a Kolmogorov-like turbulence 
spectrum is strongly disfavoured (x 2 > 100). For intermediate values of 5, good models are 
obtained for the full range in L (as in Fig. 1). For low values of 5, models with a small halo size 
L are excluded; in particular for 5 < 0.45, there is no good model with L < 15 kpc. Finally, for 
high values of S, models with a large halo L are excluded, and for 5 > 0.85, there is no good 
model with L > 1 kpc. As can be seen in Fig. 7, each value of 5 gives a different contour plot 
in the Kq/L - L plane. It appears that they all can be superimposed to a single curve by a 
rescaling Kq/L — > Kq/L x f(5), where / is a function of 5 only. For the contours displayed 
in Fig. 7, it takes the values /(0.46) = 0.51, /(0.5) = 0.62, /(0.6) = 1, /(0.7) = 1.54 and 
/(0.85) = 2.78. 

In the V a /-\/K~o - V c plane, the values of V c are shifted downward as 5 is decreased but the 
allowed range of V a / \JKq does not significantly move. Nevertheless, the allowed values for V c 
never reach 0, so that no-wind models can be excluded. 

6. Discussion and conclusion 

We obtain good quantitative constraints on the diffusion parameters from B/C data. In 
particular, there is a very strong correlation between L, Kq/L, V c , V a /\/KQ and 5. For 5 = 0.6, 
we find that 8.5 km s" 1 < V c < 12 km s _1 and 410 < V a / \J Kq < 530 (where V a is expressed 
in km s _1 and K in kpc 2 Myr -1 ). Furthermore, we show that the power law index for the 
diffusion coefficient is restricted to the interval [0.45,0.85], the best x 2 being 25.5 for 5 = 0.70. 
For any 5 in this interval, the good parameters in the Kq /L — L and V a / \JKq - V c planes can be 
straightforwardly deduced from the corresponding values for 5 = 0.6 by a simple scaling law. We 
exclude any model without a convective velocity or without reacceleration for any combination 
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of the three other diffusion parameters. 

Our conclusions could get more stringent by new measurements in the whole energy spec- 
trum for all nuclei. We emphasize that all our results were obtained using the best data, which 
are rather scarce and more than 20 year-old; new data are thus strongly needed. The AMS 
experiment on board the International Space Station will have in principle the ability to provide 
some of these data. 

The next steps of this analysis will be to study radioactive species. In particular, we expect 
the recent smili data (giving the spectral distribution of Be isotopes over a large energy range) 
to provide a new insight on cosmic ray propagation, and thus to constraint further the diffusion 
parameters. We are also investigating the standard antiproton signal, using the results described 
in this paper. 
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A. Solution for the diffusion model 

Sources and interactions with matter are confined to the thin disk and diffusion which 
occurs throughout disc and halo with the same strength is independent of space coordinates. 
Sources follow a universal form Q(E)q(r) with an associated normalization q^ for each nucleus 
(see § 3.5.2). Furthermore, we consider here only constant wind V c in the z direction so that 
adiabatic losses (third term in (1)) vanish. 



A.l. High energy limit 

The dependence in energy E is implicit, and if we add radioactive contributions localized 
in both disc and halo 

" K ad ) N j (r, z) + Y, T rL Nk + ^S(z) (qlQ(E)q(r) + ^f^N k (r, 0) - f W(r, 0)] = 0(A1) 
k=l V k=l / 

d ( d 2 19, 9 



In our geometry a solution can be found in term of Bessel functions of zeroth order taking advan- 
tage of their properties (Jackson 1975). One can expand all the quantities over the orthogonal 
set of Bessel functions { Jo(Qx)} t=1 

oo 

r 



Ni(r,z) = YNi(z)MC l -) (A3) 



i=i 



W) = 2^QiMCi-B), Qi = „ 2 j 2 ,f \ 71 -—; » P = r R A4) 

R TrR?Jt((i) J pq(p)dp 

Such an expression respects automatically one of the boundary condition N(r = R,z) = 0. We 
will use a compact notation to describe the most general form for a source term 

m k >mj 

Qj=qiQ(E)q l+ £ f^'iVf(O) (A5) 

k 



Stable progenitors N k As equations have to be valid order by order (second term of (Al) 
vanishes) 

(Qj OhVJ \ 

-j r + ^ r Ni(0)U(z) (A6) 

A fast procedure to solve (A6) is given by the following three steps: 



dz 2 



Kdz 



R 2 



+ 



rad ' 

K ' 



1. Solution in the halo (no right-hand side in eq. (A6)) with the boundary condition N-(z = 
±L) = 
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2. Integration (lim^n f^h ■ ■ ■ dz) of equation (A6) through the thin disc 6 , which gives 

2N'j(z)\ z=0 - 22V?(0)$ - 2WV/(0)§ + 0? = (A7) 

3. Put the halo solution in equation (A7) to ensure continuity beetwen the two zones. 
We finally obtain the solutions for stable progenitors in relativistic regime: 



v 7 i=o 



- . sinh 

Q? 


'si(L-z) 

2 


A j 
1 sinh 


2 





Jb(Ci^) (A8) 



= ggQ(£)& + f fc iJV* (0) (A9) 

Si = + 4|r + 4^) 1 / 2 A] = 2/if $ + F c + KSf coth(^) (AlO) 

For a primary = q J Q(E)qi, and for a pure secondary Q? = f kj N k (0). Note that 

solutions given in Webber et al. (1992) for secondary takes advantage of the primary form of 
N k (0). Since we are here interested in a shower-like (see § 3.2.2) resolution, the form given here 
is more adapted. 

(3 decay contribution from N k For all the nuclei treated here, never has more than one 
unstable contribution, so that the sum over k for iV* reduces to one term in equation (Al). 
Resolution is complicated by the localisation of this source in the whole halo. Focalising on 
this specific term, neglecting for a while primary source and classical spallative secondary con- 
tribution 2h5(z) r fcj N k (r, 0), one obtains (following the same procedure as described in the 
previous section 7 ) 

N** (r,z)= J2MCi^)x ™ d (All) 

R K ( a i _ ° ) sinh(%^ 



6 In terms of distribution (quoted in braces), denning <jo and <J\ as the discontinuities of th et 1 st order, 
remember that 



Imposing the continuity of the vertical cosmic ray current across the plane z = 0, we thus have 
at = lim^o [dNl{z)/dz\ + _] = -2N((0)% and a = 0. 

7 The contribution of these radioactive nuclei may be unimportant in some cases, but we should take it into 
account as it is the dominant process for some others. In the simple example of 10 Be— > 10 B, neglecting this channel 
would give an error of about 10% on the B flux, whereas considering that this term is only located in the disc 
would give an error of about 3% compared to the rigourous treatment given above. Notice finally that at fixed 
energy per nucleon, the rigidity depends on the nuclear species at stake. The diffusion coefficient K 1 of the child 
nucleus j is therefore different from its progenitor's one K k . The difference K? — K k tends to vanish for the 
heaviest nuclei. 
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where Sj and Aj have already been defined in (A10) and 
n = eyin (¥ck(A h\\ n =vL(JL U + . = ^£M M f A12l 

Ui — exp ^ 2 ^ fe kT))> a —2K k \K k Ki)^ K k Ki ■> a% ~ 2 \ K k pi 



A. 2. Full solution 

Modifying relation (Al) so as to take into account energy losses and diffusive reacceleration 
is straightforward since those processes take place only in the disc and not in the halo. Following 
step by step the three point procedure previously described, one gets the differential equation 8 

A\Nm = Q? - 2h^ [&L s (£0iV/(O) - ^ E (£)^/(0)} (A13) 

where bj oss (E) and K EE respectively correspond to energy losses (described in section 3.6.1 and 
defined by (21)) and energy diffusion (reacceleration term given by (20)). 



8 For the /9-unstable progenitor, the derivation is a bit more tedious but leads to the same final compact result. 
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Fig. 1. — The x 2 (adjustment to B/C data, see text) has been computed in all the parameter 
space for 5 = 0.6 (defined by K = KqTZ s ). A best x 2 is obtained for each L and Kq/L. Left 
figure displays the x 2 values in the Kq/L - L plane. It is truncated for x 2 > 50. Right figure 
shows the iso-x 2 lines of the previous surface. 




Fig. 2. — As in Fig. 1, a best x 2 is obtained for each V c and V a j \JK§. Left and right figures are 
similar to those in Fig. 1. 
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Fig. 3. — This curve displays the computed ratio of ( 10 B+ 11 B)/( 12 C+ 13 C+ 14 C) for a configu- 
ration giving a reduced xf ~ The experimental points are from HEAO-3 (solid circles), ISEE 
(triangles), imp-8 (empty circle), voyager (square) and balloons (crosses). 
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Fig. 4. — This curve displays the computed ratio of (Sc+Ti+V)/Fe for the same configuration 
as in Fig. 3. The experimental points are from heao-3 (solid circles), isee (triangles), voyager 
(square) and balloons (crosses). 
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Fig. 5. — Theses curves display the computed flux of Oxygen for the same configuration as in 
figure 3. The solid line corresponds to 6 + aoxygen = 2.68, the dotted line to 5 + aoxygen = 2.80 
(see text for details). The experimental points are from HEAO-3. 
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Fig. 6. — The parameter space has been explored for a fixed value of the halo thickness L = 3 kpc, 
varying 5 (spectral index of K) as well as the other parameters. A best x 2 is obtained for each 
S and Kq/L. This figure shows the contour lines for fixed values of x 2 . 
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Fig. 7. — Models with different values of 5 are shown. As in the previous figures, for each value 
of L and Kq/L, only the best x 2 value is retained when the other parameters V c and V a /y/Ko 
are varied. The figure in the left panel displays the contour levels for \ 2 < 40 for the indicated 
values of 5. It is possible to scale the Kq/L values by a function f(5) to superimpose the contours 
corresponding to different values of 5 (see text). This is displayed in the right panel. 
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Fig. 8. — Models with different values of 8, the coefficient diffusion spectral index, are shown. 
For each value of V c and V a /y/Ko , only the best x 2 value is retained when the other parameters 
L and Kq/L are varied. The figure displays the contour levels for x 2 < 40 for the indicated 
values of 8. 



